function [] = display_multichannel_data(imgVolume, mrprot, ifft_flag, recon_mag_flag, recon_phase_flag)
%DISPLAY_MULTICHANNEL_DATA Displays multi-channel raw MRI data loaded with rdMeas and reordered with reorder_raw_siemens_mridata.m 
%   Input : display_multichannel_data(imgVolume, mrprot, ifft_flag, recon_mag_flag, recon_phase_flag)
%   Output: No Output - displays images from reorded raw siemens MRI data.
%
%
%   recon_phase_flag: If "1", displays the phase data using sum of squares for multi-channel recon.
%   recon_mag_flag: If "1", displays the magnitude data using sum of squares for multi-channel recon.
%   ifft_flag: If "1", inverse fft the raw data before displaying it.
%
%   ***************ATTENTION****************
%   **Might not work as intended with iPat**
%   ****************************************
%
%   Author: Mahtieu Boudreau
%   Date: March 27th 2013
%   Email: mathieu.boudreau2@mail.mcgill.ca
%
%   Last updated: Devember 2nd 2013
%   Modifications: Fixed a bug in the sum of squares.


if nargin < 4
    
    recon_phase_flag = 0;
end

if nargin < 3
    
    recon_mag_flag = 0;
end

if nargin < 2
    
    ifft_flag = 0;
    disp('Data display will not be inverse Fourier transformed');
end

if nargin == 0
    
    disp('Someone forgot to type some inputs in display_multichannel_data.')
    fprintf('(???)  (???)  (???)  (???)  (???)  (???)\n') 
    return
    
end

nargin

%   Verify number of channels, set number of rows and columns for
%   concatinated image of all channels

numChannels = size(imgVolume,length(size(imgVolume)));

if numChannels == 32
    nRows = 4;
    nCols = 8;
    disp('Displaying 32 channel data');
elseif  numChannels == 8
    nRows = 2;
    nCols = 4;
    disp('Displaying 8 channel data');
elseif  numChannels == 4
    nRows = 2;
    nCols = 2;
    disp('Displaying 4 channel data');
elseif  numChannels == 1
    nRows = 1;
    nCols = 1;
    disp('Displaying single channel data'); 
else
    disp('Error:  Code not written to handle number of channels that are not equal 32, 8, 4 or 1.')
    disp('Complain about original code author for 5 minutes, then modify/fix the code yourself.')
    return
end

% Get number of slices
if length(size(imgVolume)) == 4
    numSlices = size(imgVolume,length(size(imgVolume))-1);
elseif length(size(imgVolume)) == 3
    numSlices = 1;
else
    disp('Error, can only handle input files that have 4 dimensions (3D) or 3 dimensions (single slice).')
end

% If demanded, FFT the data before displaying it
if ifft_flag == 1 
    
     % 3D or Multi-Slice
     if length(size(imgVolume)) == 4
         
         % 3D 
         if(strcmp(mrprot.tMRAcquisitionType,'3D'))
             for chIndex = 1:numChannels
                 imgVolume(:,:,:,chIndex) = ifftshift(ifftn(imgVolume(:,:,:,chIndex)));
             end
          
         % Multi-Slice
         elseif(strcmp(mrprot.tMRAcquisitionType,'2D'))
             for chIndex = 1:numChannels
                 for slIndex = 1:numSlices
                    imgVolume(:,:,slIndex,chIndex) = ifftshift(ifftn(imgVolume(:,:,slIndex,chIndex)));
                 end
             end
         end
         
     % Single slice    
     else
         for chIndex = 1:numChannels
             imgVolume(:,:,chIndex) = ifftshift(ifftn(imgVolume(:,:,chIndex)));
         end
     end
end

% Display data loop
imgCell = cell(nRows,nCols);
for iSlice = 1:numSlices
    
    % Get current slice
    if length(size(imgVolume)) == 4
        imageSlice = squeeze(imgVolume(:,:,iSlice,:));
    else
        imageSlice = imgVolume;
    end
    
    % Useful Trick - input images in cell, then use cell2mat to concatinate the
    % images in a single 2D matrix.
    for iImage = 1:nRows*nCols

       imgCell{iImage} = imageSlice(:,:,iImage);

    end

    bigImage = cell2mat(imgCell);

    % Plot concatinated images of all channels for this slice.
    figure(1)
    imagesc(abs(bigImage));
    title(sprintf('Slice %i', iSlice));
    truesize

    % Display magnitude and/or phase of sum of squared reconstructed data. 
    if (recon_mag_flag == 1) || (recon_phase_flag == 1)

        for ii = 1:numChannels
            if ii == 1
                sumsqrImage=abs(imageSlice(:,:,ii)).^2; 
            else
                sumsqrImage=sumsqrImage+abs(imageSlice(:,:,ii)).^2; % Sum of squares of the *complex* data.  
            end
        end
        
        % Display magnitude
        if recon_mag_flag == 1
            figure(2)
            title(sprintf('Magnitude image of slice %i',iSlice))
            imagesc(abs(sqrt(sumsqrImage)))
            truesize
        end
        
        % Display phase
        if recon_phase_flag == 1

            figure(3)
            title(sprintf('Phase image of slice %i',iSlice))
            imagesc(angle(sqrt(sumsqrImage)))
            truesize
        end
    end
    
    
    drawnow;
    pause(2);
end



end

